clear;clc;
global theta e12 eg ey;

Nblack = 607;
Nwhite = 3973;

rng(10010);
N_sim=100000;

wparam=[0.174027339165297;0.187646010085453;1.52703723856717;0.578390833256311;1.47397496404516;0.472254474029496;0.539039975907426;0.508222124622837;0.557092030919351;1.68356685222368;-0.456906445371444;1.53322800556343;0.520506916649071;0.554621783550490;0.495990534196635;0.550405334264551;0.503398932635956;0.149831038482197;1.21750922668894;0.774068381366255;-0.0943955231130527;-0.454412676623209;0.225964287920840;-7.47819116937451e-05;0.226994852644431;0.162646835240202];
bparam=[-0.752900651375512;-0.625236744850452;0.955164978663548;0.470520928662698;0.940898873680747;0.457615100857030;0.522441225000711;0.796501681223757;0.526383943012879;1.38074364109675;-0.825533240871413;1.01277501698395;0.503940752081312;0.230058029470066;0.270326488861762;0.435689411502657;0.141913310661526;-0.625025361512185;0.665429321866805;0.834583035216824;-0.257622993708532;-0.213896500811884;-0.532316091400443;-0.518906855782177;0.0476680272326578;0.309398571910553];


%% Distribution of \theta
  
theta = normrnd(0,1,N_sim,1);
w.theta = theta*wparam(8);
b.theta = theta*bparam(8);
w.eg = normrnd(0,1,N_sim,1)*wparam(20);
b.eg = normrnd(0,1,N_sim,1)*bparam(20);
ey = normrnd(0,1,N_sim,1);


%% DIFFERENT CASES

%% Both White Teachers
param = bparam; theta = b.theta; eg=b.eg; e12 = normrnd(0,1,N_sim,2);
a=simulation(param,1);
[f_bo,x_bo]=ksdensity(a(:,9),'function','cdf','bandwidth',0.001);

% Whites, SR
param = wparam; theta = w.theta; eg=w.eg;
a=simulation(param,0);
[f_w,x_w]=ksdensity(a(:,9),'function','cdf','bandwidth',0.001);

% due to cy, theta and GPA
param = bparam; theta = w.theta; eg=b.eg; e12 = normrnd(0,1,N_sim,2);
param(11)=wparam(11);
param(4)=bparam(4)-bparam(11)+wparam(11);
param(9)=bparam(9)-bparam(11)+wparam(11);
param(18:20)=wparam(18:20);
a=simulation(param,1);
[f_bwcytg,x_bwcytg]=ksdensity(a(:,9),'function','cdf','bandwidth',0.001);

% Outcome due to teacher expectations
param = bparam; theta = w.theta; eg=b.eg; e12 = normrnd(0,1,N_sim,2);
param(4:5)=wparam(4:5);
param(16)=wparam(16);
param(9:10)=wparam(9:10);
param(4)=wparam(4)-wparam(11)+bparam(11);
param(9)=wparam(9)-wparam(11)+bparam(11);
param(17)=wparam(17);
% cy param 
param(11)=wparam(11);
param(4)=bparam(4)-bparam(11)+wparam(11);
param(9)=bparam(9)-bparam(11)+wparam(11);
param(18:20)=wparam(18:20);
a=simulation(param,0);
[f_bwte,x_bwte]=ksdensity(a(:,9),'function','cdf','bandwidth',0.001);


% Give BLACKS the WHITE PRODUCTION FUNCTION VALUES
param = bparam; theta = b.theta; eg=b.eg; e12 = normrnd(0,1,N_sim,2);
param(13:15)=wparam(13:15);
a=simulation(param,1);
[f_bwy,x_bwy]=ksdensity(a(:,9),'function','cdf','bandwidth',0.001);



figure
subplot(2,2,1)
plot(x_bo,f_bo,'k-.',x_w,f_w,'k');
legend('Blacks','Whites','location','southeast');

xlabel('Pr(Y=1)')
ylabel('CDF')
xlim([0,1]);ylim([0,1]);

subplot(2,2,2)
plot(x_bwcytg,f_bwcytg,'k.-',x_bo,f_bo,'k-.',x_w,f_w,'k');
legend('Blacks, White \theta, G','location','southeast')

xlabel('Pr(Y=1)')
ylabel('CDF')
xlim([0,1]);ylim([0,1]);


subplot(2,2,3)
plot(x_bwte,f_bwte,'k--',x_bo,f_bo,'k-.',x_w,f_w,'k');
legend('Blacks, White \theta, G, T_j','location','southeast')

xlabel('Pr(Y=1)')
ylabel('CDF')
xlim([0,1]);ylim([0,1]);


subplot(2,2,4)
plot(x_bwcytg,f_bwcytg,'k.-',x_bwte,f_bwte,'k--',x_bo,f_bo,'k-.',x_w,f_w,'k');
legend('Blacks, White \theta and G','Blacks,  White \theta, G, T_j.','location','southeast')

xlabel('Pr(Y=1)')
ylabel('CDF')
xlim([0,1]);ylim([0,1]);

saveas(gcf,'white_teacher.jpg')





figure
plot(x_bo,f_bo,'k-.',x_w,f_w,'k');
legend('Blacks','Whites','location','southeast');
xlabel('Pr(Y=1)')
ylabel('CDF')
xlim([0,1]);ylim([0,1]);
saveas(gcf,'cf_1.jpg')

figure
plot(x_bwcytg,f_bwcytg,'k.-',x_bo,f_bo,'k-.',x_w,f_w,'k');
legend('Blacks, White \theta, G','location','southeast')

xlabel('Pr(Y=1)')
ylabel('CDF')
xlim([0,1]);ylim([0,1]);
saveas(gcf,'cf_2.jpg')

figure

plot(x_bwte,f_bwte,'k--',x_bo,f_bo,'k-.',x_w,f_w,'k');
legend('Blacks, White \theta, G, T_j','location','southeast')

xlabel('Pr(Y=1)')
ylabel('CDF')
xlim([0,1]);ylim([0,1]);
saveas(gcf,'cf_3.jpg')

figure



plot(x_bwcytg,f_bwcytg,'k.-',x_bwte,f_bwte,'k--',x_bo,f_bo,'k-.',x_w,f_w,'k');
legend('Blacks, White \theta and G','Blacks,  White \theta, G, T_j.','location','southeast')

xlabel('Pr(Y=1)')
ylabel('CDF')
xlim([0,1]);ylim([0,1]);

saveas(gcf,'cf_4.jpg')
clear;
clc;